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Abstract 

This paper addresses the deconvolution of an image that has been obtained by superimposing 
many copies of an underlying unknown image of interest. The superposition is assumed to not 
be exact due to noise, and is described using an error distribution in position, orientation, or 
both. We assume that a good estimate of the error distribution is known. The most natural 
approach to take for the purely translational case is to apply the Fourier transform and use the 
classical convolution theorem together with a Weiner filter to invert. In the case of purely rotational 
deblurring, the similar Fourier analysis is applied. That is, for an blurred image function defined on 
polar coordinates, the Fourier series and the convolution theorem for the series can be applied. In 
the case of combinations of translational and rotational errors, the motion-group Fourier transform 
is used. In addition, for the three cases we present the alternative method using Hermite and 
Laguerre-Fourier expansion, which has a special property in Fourier transform. The problem that 
is solved here is motivated by one of the steps in the cryo-electron-tomographic reconstruction of 
biomolecular complexes such as viruses and ion channels. 
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1 Introduction 

In single particle Cryo- Electron-Microscopy (Cryo-EM), the goal is to reconstruct the 3D shape of 
large biomolecular complexes from projection data. In particular, many essentially identical copies 
of a complex of interest are embedded in a thin layer of vitreous ice at randomized (and unknown) 
orientations. If we consider this thin film of ice to be in the x-y plane in the lab frame, then an 
electron beam takes projections of the density of the embedded biomolecular complexes along the 
z direction. The goal is to reconstruct the three-dimensional density of the complex (which defines 
its shape) from these projection images. The difference between this problem and medical image 
reconstruction is that the projection directions are unknown a priori. 

The signal-to-noise ratio in such measurements can be quite high [6]. Therefore projections 
corresponding to the same (or quite similar) projection directions are grouped together and su- 
perimposed. In doing so, the random noise of the background has a tendency to cancel, and the 
features of interest in the projections reinforce each other as the number of superimposed projec- 
tions becomes large [6] . This averaging technique may not be used in all types of 3D reconstruction 
of the electron microscopy. However, it is still useful for the analysis of classes of 2D projected 
images, the raw data of which contains a large amount of noise. 

One problem is that the superposition of images might not be exact. This results in a blurring 
relative to the true underlying image of interest. To get a sense of this, let us consider the following 
image model including zero- mean ergodic white noise [6] . 

p{x,t) = po(fft"^x) +n(x,i). (1) 
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Here gt is the homogeneous transformation in SE(2) ( the Lie group describing the translational and 
rotational motion), n is the noise, x e is the planar position of points in each image, and t e R"*" 
is an artificial time variable used to order the images. If there is no noise term, it is intuitively 
clear that the appropriate matching of two images, p(x, ij) and p{x,tj) occurs at gt^ — gt. and the 
superposition is (p(x, ti) + p(x, tj))/2 — /o(g'i~^x). However, if we have the noise term, the matching 
of many data images produce various gt- The superposition of them will be in the form of 

-^p(x,t,)~^^Po(5t:'x)= / {-Y^8{g-^og)\p^{g-^^). 
i=l i=l ''^ \ i=l / 

The first equality assumes that the noise term is approximately canceled out during the supcrim- 
position, and the second equality shows that this superposition is a convolution on the group of 
rigid-body motions of the plane, G — SE(2). The right hand side is essentially the blurred version 
of po(x) depending on the distribution of gt^. 

We therefore seek to solve the following inverse problem: Given a blurred image, 7(x), that de- 
scribes the optimal superimposition of many experimentally-obtained projection images, and given 
an estimate of the probability density function describing the error distribution in the alignment of 
these superimposed images, f{g), we seek to find the deblurred image p(x). This is expressed as 
the solution to the problem: 

/ f(g)p(g-'-^)dg^^(^). (2) 
Jg 

Here G is the group of transformations involved in alignment, g ■ x denotes the group action of 
G on R^, and dg is the associated invariant integration measure for that group [5]. In this paper 
we consider three cases: (1) G = (R^,-|-), the translation group in the plane; (2) G — SO{2), 
the rotation group in the plane; and (3) G = SE(2), the Euclidean motion group of the plane. 
Explicitly, in these three cases we have 
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/oo poo 
/ fi{yi,y2)p{xi-yi,X2-y2)dyidy2 = -fi{xi,X2), (3) 
■oo J —oo 

j-2-K 

/ /2(^)p(a;i cos^ + X2sm^; — xi sin^ + ^2 cos^)d^ = 72(xi, X2), (4) 
Jo 



/oo /"OO Z'.^TT" 
/ / /3(yi, Z/2, - Vi) cos 6* + (a;2 - ^2) sin 9; -{xi - yi) sin 6* + (^2 - ^2) cos e)dyidy2de 

00 J —00 J 

(5) 

= 73(a;i,a;2). 

As a model for the functions fi, we will assume appropriate concepts of Gaussian distributions. 
Recall that the diffusion equation on the line, 

dt dx'^ ' 

subject to the initial conditions f{x, 0) = S{x), has the solution 

= ^^e-^'/^*. (6) 



The solution to the uniform planar diffusion equation: 

dt dx\ dx\ ' 

can be written as 

X2, t) = /(xi, t)J{x2, t). 
Likewise, the diffusion on the circle can be viewed as a folded normal distribution: 



f2{e,t)^ f{e-2nn,t). 



It is useful to note that this can be expressed alternatively as a Fourier series: 




A model for combined translational and rotational error is then 



f3{Xi,X2,0;ti,t2) = fl{Xi,X2,ti)f2{9,t2) 



where the small values of ti and t2 can be chosen separately to describe different amounts of 
translational and rotational error, as well as to account for the fact that the units of measurement 
are different for translations and rotations. 

2 Deconvolution of Motion- Averaged Images Using Fourier 
Transform 

In this section we address how to solve each of the three deconvolution problems in ([3]) - (|5]) using 
Fourier transform. 

2.1 The 2-D Fourier Transform for Translational Deconvolution 

The natural tool to use to solve the deconvolution problem in ^ is the Fourier transform. The 
Fourier transform in two dimensions is written in Cartesian coordinates as 




and the inversion formula is 




Fourier transform of the distribution function, fi{xi, X2,t) in (h| is fi{uj) = e ('^i+'^z)*. 
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The convolution theorem then converts dSl) to a problem in Fourier space of the form 



fi{uj)p{uj) = 71 (^), 



which is inverted after regularization as 



pM = 7iH/iH/(6 + |/iMr). 



(7) 



The regularization parameter, e, is a very small positive number that is introduced to handle zeros 
of the Fourier transform. In fact, for the Gaussian distribution of interest in our problem, there 
are no zeros, but in both real and Fourier space the tails of the distribution can approach zero at 
points sufficiently far from the origin. ([T]) is nothing more than the well-known Wiener filter. 

2.2 Deconvolution of Purely Rotational Misalignment 

If the image functions are defined on polar coordinates, Q can be rewritten as 



^0 

where xi = rcos0 and X2 = rsin0. If we fixed the value of r, Q becomes the convolution of the 
two functions on a circle as 

/ Me)p^'H<P-o)de = ^!{\<j>), (8) 

Jo 

where p^'^\4> — 9) = p{r, 4> — 9) and 72'''' (0) = 72(^^5 0)- 

The Fourier series expansion of a function defined on a circle gives 




^ 00 



—00 



where 



n 
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The Fourier transform of the distribution function, f2{6,t) in (|4|) is (/2)n = *. 

The convolution theorem of Fourier series converts ([s]) to the problem in Fourier space of the 
form 

(/2)n(p(^))n = 

As in the case of the translational deconvolution, the inversion with regularization is 



=(¥^))n(/2)n/(6+|(/2)nP). (9) 

2.3 Deconvolution of Combined Translational and Rotational Blurring 

In order to solve the full motional deconvolution problem, the appropriate concept of Fourier trans- 
form is required. In particular, since /s is a function on the group of rigid-body motions of the 
plane, SE(2), and a function on can be viewed as a function on SE(2) that is constant over the 
orientational variable, ([s]) can be viewed as a convolution on SE(2). We therefore review here the 
group SE(2) and the associated Fourier analysis. 

2.3.1 Representation Theory of The Euclidean Motion Group of the Plane 

Each element of SE(2) is parameterized in either rectangular or polar coordinates as: 

(cos 9 — sin 9 ai 
sin 9 cos 9 a2 
1 

or 

(cos 9 — sin 9 a cos (p 
sin 9 cos 9 a sin (p 
1 

where a = 11 all. 
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A irreducible unitary representations of SE(2) (see [5,10,11] for general definition) can be viewed 
as infinite dimensional matrices, U{g,p) with elements expressed as 



Ur 



where Jy{x) is the z/*'* order Bessel function and m and n range over all integer values. 

From this expression, and the fact that U{g,p) is a unitary representation, we have that: 

u„,n{g~^{a,(f),e),p) =u;;^^{g{a,(j),e),p) = 

Ur^m{9iaA,0),p) = 2"-™e^["^^+("-"^)'^l J„_„(pa). (11) 

These matrix elements are related by the symmetries: 



Umn{g,P) = (-1)" "'u^m,-n{9,P), (12) 

Umn{g{-a',4>,0),p) = Umnigia,(j)±n,e),p) = (- 1 )"""«„,„ (^((a, 0, 61) , p) (13) 

and 



(-l)'""'^n„,„((7(a, <P-e, -e),p) = Unmigia, 0, 0),p). (14) 



The equality in (14) follows from (11) and (13). 



2.3.2 The Fourier Transform for the Euclidean Motion Group of the Plane 

The Fourier transform of a sufficiently well-behaved function on SE(2), and the corresponding 
inverse transform are defined as: 

Hf) = hp)= I f{9)U{g-\p)dg 
Jg 

and 



/• oo 

:^~\f) = f{9)= / trace 
Jo 



{f{p)U{g,p))pdp. 



As with the Fourier transform of functions on 



A proof that these identities hold is given in [10]. The fact that the inverse transform works depends 
on {U{g,p)} being a complete set of irreducible representations, and the fact that it is unitary allows 
us to write U{g~^,p) = instead of computing the inverse of an infinite dimensional matrix. 

The matrix elements of the transform can be calculated using the matrix elements of U{g,p) 



defined in (10) as: 

fmn{p)= / f{g)Umn{g'\p)dg. (15) 

Jg 

Likewise, the inverse transform can be written in terms of elements as: 



{g,p)pdp. 



2.3.3 Regularized Deconvolution of Motional Deblurring in SE(2) 

Given motional blurring expressed in ([2]) when G = SE(2), the result becomes Q. This IS a 
convolution on the motion group. The result can be solved by applying the motion group Fourier 
transform to yield: 

p{p)Mp) = 73 (P)- (16) 
The functions p{p) and 73 (p) are row vectors. The direct matrix inversion would be 

Kp) = %{p)[kp)]~'- (17) 

However, if the matrix f^ip) becomes singular, then this needs to be regularized. The procedure 
for doing this is explained in [4] , and involves the computation of a weighted least-squares pseudo- 
inverse. 



In the current context, we can compute the entries of fsip) analytically. For the time being, we 
drop the subscript '3', and write the function in polar coordinates as 

^ oo 

/(r,<^,^;ti,t,) = -4-e-^V4*. ^ ^-k^t,^^ke_ 

oTT Zi 

fe=— oo 

We use the fact that in polar coordinates dg = rdrd(j)d9 and the above function is independent of 
(f) (which will result in a diagonal SE(2) Fourier transform matrix). Computing the SE(2) Fourier 
transform of the distribution, we find: 

fmn{p) = ^5mn (j^ '^'^ Mpr)rdr\ e'"^"'^ = 5^„e-P'*^e-"^'*^ (18) 



Therefore, the matrix, fz{p) in (17) is diagonal. Its inversion is the simple inversion of scalar values 



and the inversion with the regularization parameter is the same as that in the previous two cases. 

3 Deconvolution of Motion- Averaged Images Using Her- 
mite and Laguerre functions 

Even though the Fourier transform is a good way to solve the deconvolution problem, its imple- 
mentation needs several manipulations of data such as interpolation. We will discuss the details in 
the next section. In this section, we develop the alternative method of deconvolution using Hermite 
and Laguerre functions. We utilize special properties that Hermite and Laguerre-Fourier expansions 
have in the Fourier transform. In this method, (|3]) and (|5]) are solved in the Fourier space, while 
(|4]) is solved in the real space. 

Hermite function, /i„(x) is an eigenfunction for the Fourier transform as 



K{x)e-"^''dx = V27r(-z)"/2,„(u;), 
where hn{x) is the Hermite function defined as 

Kix) = -Hn{x)e~^"/\ 
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where s„ = •\/2"n!y^ and Hn{x) is Hermite polynomial, which is generated by the Rodrigues 
formula 



2 (i" ^ _^2, 



Hn{x) = i-lTe^ — (e 

This property gives straightforward analytic solution, when the image is defined as Hermite expan- 
sion and the Fourier method in the previous section is applied. 

While Hermite expansion is good for a function defined on Cartesian coordinates, Laguerre- 
Fourier expansion looks better for a function defined on polar coordinates. Mathematically, the two 
expansions can be converted to each other [1]. 

The associated Laguerre polynomials arc generated by the Rodrigues formula. 

The associated Laguerre polynomials are orthogonal over [0, oo) with respect to the weighting 
function x^e~^, 

f k —X r k / \ T k ( \J ~^ ^) ■ X 

Jo 

Using the Laguerre polynomials and Fourier basis, we can define the basis function on two-dimensional 
polar coordinate as follows [1] [2] . 



e '""^ 



where (m — |n|) and (m -|- |n|) are even numbers. For convenience, we sometimes divide it into two 
parts as 

Xm,n(r,(/)) =ym,n(r)^„(0), 

where 



[(m — 


n\)/2]\ 


7r[(m -|- 


\n\)/2]\ 
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In this section we will show how to solve the aforementioned deconvolution problem using the 
Hermite and Laguerre functions. The appropriate expansion will be chosen for the three cases. 

3.1 The Translational Deconvolution using Hermite expansion 

Hermite expansion of an image function is 

oo oo 

p(Xi,X2) = ^^Pmnhm{Xl)hn{x2), 
m=0 n=0 

where 

Pmn= / p{Xi,X2)hm{Xi)hn{x2)dXidX2. 

When an image function can be expressed as a truncated Hermite expansion with large N, the 
function is written as 

N N-m 

p{xi,X2) = XI 2Z Pmnhm{Xl)hn{x2). (19) 
m=0 n=0 

The Fourier transform of the image function is 

N N-m 

p{uJi,UJ2) = ^Y1 Pmn{'2n){-i)'"+''hm{uJl)hn{uJ2). 
m=0 n=0 

Since the Fourier transform of fi{xi,X2,t) is f 1(001,002) — e^^'^i+'^a^*, the convolution theorem gives 

N N-m 

7i(a;i,u;2) = ^ E Pmn(27r)(-i)"^+"/i^(u;i)/i„(a;2)e-('^'+^^')*. (20) 

m=0 n=0 

On the other hand, since H^{x) is an m'th order polynomial, it can be rewritten as 



a 



k=0 



where am,k{ci ^) is an appropriate coefficient relating Hermite polynomial and its scaled version. 
Using this expression, we can have 

m 

hm{uj)e~'^ * = J^am,fc(a"^) — hk{auj), 
k=o 
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where a = \/2t + 1. Therefore, 

A' N—m m n 

7i(t^i,t^2) = X] Pmn(2vr)(-i)™+"^^am,fc(a"^)a„,z(a~^) hk{auji)hi{auj2). 

m=0 n=0 fc=0 1=0 

We can reorder the summations and have 

N N~k / N-l N~m 



li[yJi,^2} = X^X] X] Pmn(2vr)(-i)™+"a„^fc(a ^)a„,i(« ^) — — hk{auji)hi{auj2). (21) 



_u , _u Sk Si 

IJmnK^'^ — ^m,ky^ 

k=0 1=0 \m=k n=l 

Its inverse Fourier transform is 

^k+l 



N N-k /N-l N-m 

7i(xi,X2) = X^X^ Pm„(27r)(-i)"'+"a„,fc(a"^)a„,z(a"^) — — j ^ — -hk{xi/a)hi{x2/a). 



k=0 1=0 \m=k n=l ™ " 

Conversely, if we can have a truncated Hermite expansion for a blurred image as 

N N~k 

7i(xi,X2) = ^^^kihkixi/a)hi{x2/a), (22) 

k=0 1=0 



then its Fourier transform is 



N N-k 27ra^ 
71 (^1' ^2) = 5^ XI 

'Iki .^_|_; hk{0'UJi)hi[auj2) ■ (23) 

k=0 1=0 



Equating (20) and (23) on a various samples on {cu^p\uj^'^^) gives 

{EHU)R{EHUf = H^x G X H^, 

where Em,n = ^m.ne"*^'^''"')', i?m,„ = /l„_i(u; ("")), {Ha)m,n = hn-liau'^'^^), Um,n = 5m,n(-0'"5 Rm,n = 

Pm-i,n-i and Gm,n = 7m-i,n-i- In Order to get R, which is the Hermite coefficients for the deblurred 
image, we should examine the inversion of the matrices. 

Inversion of U is given by = 6m,ni"^- Pseudo-inverse of H is given by = {H^H)~^H^ if 
the sampling points, o;^^-* are chosen appropriately as shown in the previous work [1]. Inverting E 
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needs regularization because inverse of e"*^'^''"'^)^ may be unstable with large value of uj^"^\ Therefore, 

^m,n — '^m,n(l/(6~**''^*'"'''^ + ^) with a Small number, e. 
Now we have 

R = U-^H+E+Ha xGx {U'^H+E+Haf. 
3.2 The Rotational Deconvolution using Laguerre-Fourier expansion 

A 2D function defined on polar coordinate can be expressed as 

oo m 
m=0 n=—m 

where 

Pmn= / p{r,(f))xmn{r,(f))rdrd(f) 
If the function can be expressed as a truncated Laguerre-Fourier expansion with large A^, we have 

N m 

p{r,<P) = Y,Y. PmnX*mnirA)- (24) 

m=0 n=—m 

Note that the integer variable n increases by multiples of 2. 

If the image function in (|4]) is defined on polar coordinates, the convolution can be rewritten as 

h{d)p{r,<P-d)dd = ^2{rA), 

with the coordinate conversion, xi = rcos9 and X2 = rsinO. If we use the truncated Laguerre- 
Fourier expansion for p, we have 



„2n -,00 N m 

' fc=— 00 m=On=—m 



The left hand side can be computed as 

00 N m p2tt 00 N m 

2i~ 



^ 00 N m „2tt 00 N m 

^ E E E ^''^'PmnyUr) / e^^'e'^^'^-'^de = E E E (^''''(^''"^Pmnymn{r)Sk, 



k=—oo m=0 n=—m k=—oo m=0 n=—m 

14 



N m N m 

m=0 n=—m m=0 n=—m 

Therefore the blurred image is 

A'' m 

<P), (25) 

m=0 n=—m 

where 

Tmn — Pmn^ 

This means that the convolved (blurred) image of a truncated Laguerre- Fourier expansion with 
purely rotational motion is also a truncated Laguerre-Fourier expansion with the same truncation 
hmit. The only difference is that the coefficients are scaled. Once we compute the Laguerre-Fourier 
coefficients(7^„) of the blurred image, the Laguerre-Fourier coefficients of the deblurred image is 
given by 

P^n = 7mn{V(e""'* + e)} (26) 

with regularization. 

3.3 The Translational and Rotational Deconvolution using Laguerre- 
Fourier expansion 

In this section, we develop a method for deconvolution of the blurred image with the translational 
and rotational motions. We utilize the special property of the Laguerre-Fourier expansion in SE(2) 
Fourier transform. 
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3.3.1 Fourier transform of the Laguerre-Fourier expansion 

When a function on polar coordinates, p{r, 4>) is defined as a truncated Laguerre-Fourier expansion 
as 

N k 

p{rA) = ^Y^ PkiXliir,(p), 

k=0 l=-k 

its Fourier transform in SE(2) is 

r"27r /'2tt poo 



pZn pZn roo 

Je=o J (t>=0 Jr=0 

^ ^ /"27r /"27r poo 

fc=o i=-k -^^=0 -^'^=0 '^'•=0 

^ ^ foo i'2'K /•2tt 

= i"— Y.Y.P'^n yki{r)Jm-n{pr)rdr / e^'-^e'^^-^^^^^ / e'^^'dO 

k / poo \ 

= 47r^i"""* X] XI '^^^^ ( / yki{'r)Jm-n{w)'rdr ) (5, _„(5„,o 

N k , „oo \ 

= 47r^i" X X '^^M / yki{r)Jiipr)rdr ) 5; _„(5^,o 
ifc=0 /=-fe ^-^^=0 / 

On the other hand, we have a useful identity in [3] and [5] as 

poo 

/ {ar)"'L';^{a^r^)e-'^'''/^Jm{kr)rdr = (-l)"Q;-'(A;/Q;)'"L;^(A;7Q;')e-'='/'"' 
Jo 

Using this identity, we can have 



r=0 



yki{r)Ji{pr)rdr^{-lY''-\^\y^J 



7rP + |Z|)/2]! 
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Therefore, 



N k 



Pmn(p) = 47rV^ ^ Pm(-1) 2 ?/fc,i(p)(5i_„(5m,0 = 47r^^" pk-n{~l) ^'^ yk~n{p)Sm,0, (27) 

fc=0 l=-k k=\n\ 

where [n/2] = n/2 if n is even and [n/2] = {n — l)/2 if n is odd. 
3.3.2 Deconvolution using the Laguerre-Fourier expansion 



In Section 3.1, we noticed that the translational blurred version of (19) is (22). This scahng effect 



on the domain appears in the polar coordinates as follows: If the original image is given as 

N m 

P(^,0) = X1 XI P'""^mn(^,0), (28) 



then its translational motion blurring is 



7(r,, 



m=0 n=—m 



N m 
m=0 n=—m 



mn \ 5 r 7 5 

a 



(29) 



since the (28) and (19) are equivalent under the simple coordinate relation, r = cos^ and r = sini 
[!]■ 



In Section 3.2, the rotational blurred version of (24) retains the structure of the truncated 



Laguerre-Fourier expansion. In other words, if the original image is expressed as a truncated 
Laguerre-Fourier expansion, then its rotational blurring gives a truncated Laguerre-Fourier expan- 
sion with the same truncation limit and domain. 



Therefore, we can conclude that the motion blurring in SE(2) of (28) has the structure of (29) 



because the full motion in SE(2) can be decomposed into translation and rotation. The details of 
the commutativity of the two motion blurring will be shown in the appendix. 



From (27) and the convolution theorem, we have 

r N — n 
I 2 

47r2i"5^_o XI Pk,-n{-l)^yk,-n{p)e~^^^^e 

k=\n\ 



7m,n(p) = \p{.P)h{.P) 
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Also, we can have the Fourier transform of (29) directly as 



7m,n{P) = 471" i Omfl ^ _„(-!) 2 a yk-n{ap). 

k=\n\ 



Equating the two expressions of 7m,n(p) gives 



k=\n\ 



^ 7fc-n(-l)2a^l/fc,n(ap) 

k=\n\ 



because ym,n = Vm-n- Since this should holds for any choice of p with a fixed n, we can have the 



following matrix expression as in Section 3.1 



WYJ X r = a^YaJ x g, 
where Wk,i = 6k,ie-(-P'''^''^e-^''\ Jk,i = 

P\n\-n P\n\+2,-n ' ' ' P2[^]+r; 
7|n|-n %\+2-n ■■■ 72[^]+r 



Y 



and 



y\n\+2,n 



J^^ = J and pseudo-inverse of F is given by Y^ = (Y'^Y)^^Y^ if the sampling points, p's are 
chosen appropriately as shown in the previous work [1], which guarantee the stable inversion of Y. 
Since inverting W needs regularization, Wj^i = (5fc,/(l/(e^*^''P*'''-'^e~*2"^ + e) with a small number, e. 
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Figure 1: Original image. The original image is from http://www.med.univ-angers.fr/discipline 
/ radiologie /Intlatlas /t laxl 1 .html 

Now we have 

r = a^JY+W+YaJ x g, 

Note that this holds for a fixed value for n. Therefore, applying it for a different n gives the full 
Pm,n, which is the Laguerre- Fourier coefficients for the deblurred image. 

4 Numerical Examples 

We generate artificially blurred data as follows; We sample the amount of motions from Gaussian 
motion distributions. Then the original image is shifted by the sampled motion. We prepare a 
set of the shifted images and average them to have our 'experimental' blurred image, 7(x). With 
knowledge of the motion distribution f{g)., we examine if inverting will give back a good estimate 
of the original. The original image and our artificially blurred images are shown in Fig. [l] and Fig. 
[2| respectively. We averaged 100 shifted images for each blurred image in Fig. |2j 
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(a) Blurred image with translational (b) Blurred image with rotational (c) Blurred image with translational 
motion motion and rotational motion 

Figure 2: Blurred image 
4.1 Case 1 : Deblurring of translational blurring 



The deconvolution method described in Section 2.1 can be implemented by discrete Fourier trans- 
form(DFT). In order to capture the detail of the motion distribution, fi, we need to sample /i on 
a very fine grid, because the distribution function is highly concentrated near the mean value. To 
match the consistency DFTs of the samples of /i and the discrete image data, we have to resam- 
ple the image on the fine grid on which we sample /i. This resampling is done by simple linear 
interpolation of the image. 

On the other hand, in order to implement the deconvolution method using the Hermite expansion 



described in Section 3.1 , we have to have an appropriate truncated Hermite expansion of the blurred 
image. The process to obtain a truncated Hermite expansion of a discrete image was developed 
in [1]. Essentially, a truncated Hermite expansion optimally fit to the image function is obtained. In 
contrast to the Fourier method, this method gives the interpolation values automatically, because 
the fitted truncated expansion is a continuous function. Furthermore, in the formulation, the blurred 
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(a) Deconvolution of the translational (b) Deconvolution of the translational 
blurred image using Fourier transform blurred image using Hermite expan- 
sion 

Figure 3: Deconvolution of the translational blurred image 
image function retains the structure of the original expansion. After having the truncated Hermite 



expansion, the estimation of the original image can be obtained by the method in Section 3.1 

Numerical results of the two methods for the deconvolution of the translational blurring are 
shown in Fig |3j 

4.2 Case 2 : Deblurring of rotational blurring 

As the implementation of the Fourier method for deconvolution of the translational motion blur 



in Section 2.1, the Fourier method deconvolution of the rotational motion blur in Section 2.2 also 
needs resampling of the blurred image on a fine and polar grid, because the distribution function, 
/2 is highly concentrated and the formulation is done on polar coordinates. Thus, interpolation of 
the blurred image on a fine polar grid is performed to match the DFTs of the distribution and the 
blurred image. 



For the method using the Laguerre-Fourier expansion in Section 3.2, firstly we have a trun 
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(a) Deconvolution of the rotational (b) Deconvolution of the rotational 
blurred image using Fourier transform blurred image using Laguerre- Fourier 

expansion 

Figure 4: Deconvolution of the rotational blurred image 

cated Laguerre-Fourier expansion for the discrete blurred image. Specifically a truncated Hermite 
expansion for the image is obtained first and then we convert it to a truncated Laguerre-Fourier 



expansion [1]. The Laguerre-Fourier coefficients of the deblurred image can be computed by (26). 

Numerical results of the two methods for the deconvolution of the rotational blurring are shown 
in Fig m 

4.3 Case 3 : Deblurring of combined translational and rotational blur- 
ring 

In order to have Fourier transform in SE(2) of the blurred image, we resample it on fine polar grid. 



With this sampling, we can compute the Fourier transform by the discretized version of (15). The 



Fourier transform of the original image can be estimated by (17). 

As in the case 2, once we have a truncated Laguerre-Fourier expansion for the blurred image, 
we can obtain the Laguerre-Fourier coefficients of the deblurred image by the method in Section 
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(a) Deconvolution of combined mo- (b) Deconvolution of combined mo- 
tional blurred image using Fourier tional blurred image using Laguerre- 
transform Fourier expansion 

Figure 5: Deconvolution of combined translational and rotational blurred image 



3.3 which is given as matrix manipulations. 

Numerical results of the two methods for the deconvolution of the translational and rotational 
blurring are shown in Fig [5j 

5 Conclusion and discussion 

In this work, we have shown our method to restore the 2D blurred images that were generated by 
the translational motion, the rotational motion, or both. In our formulation, the blurring process 
can be expressed as a convolution of an original image and the motion distribution function. Since 
the convolution can be interpreted as the multiplication in the Fourier space, the deconvolution ( 
i.e. deblurring ) is the simple inversion process in that space. 

We apphed this concept to the three cases. We used the Fourier analysis in SE(2) for the 
combined translational and rotational motion blur, which is defined using the representation theory 
of the Euclidean motion group of the plane. 
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The Fourier analysis has the strong advantage that the convolution in the original space can be 
replaced with the multiplication in the Fourier space, which is much simpler to manipulate. However 
the implementation of it needs more manipulation of data, because the image is the function defined 
only on a discrete grid domain, while the formulation is done in the continuous domain. Although 
the discrete Fourier transform is a good way to do it, resampling is required before the transform as 
described in Section |4j Shortly, the resampling is needed for matching the resolution of the image 
function and the distribution function. In addition to the resolution matching, we need another 
resampling process for the cases of the rotational motion blur and the combined motion blur, since 
the image function on the polar grid is required in these cases. 

To overcome this, we proposed an alternative method using the Hermite and Laguerre-Fourier 
expansions. Once the 2D truncated Hermite expansion optimally fitted to the given discrete image 
is obtained, we can directly utilize the results of the derivation without introducing discrete version 
of it, which is essential in the Fourier method. In this case, resampling on a finer grid comes 
naturally, because our expansion is already a function on a continuous domain. On the other hand, 
we already know that the interconversion of the Hermite and Laguerre-Fourier expansions is possible 
losslessly and it can be viewed as the coordinate conversion from Cartesian to polar coordinates [1] . 
Therefore, once we have the Hermite expansion for an image, the corresponding expression on the 
polar coordinates, which is exactly the Laguerre-Fourier expansion, can be obtained easily. Since the 
two expansions have the special property under the Fourier transform, which is that they retain their 
structure under the Fourier transform, the expansions enable the straightforward implementation 
of the deconvolution method. 

The computational complexity of the interconversion of the two coordinates is still high(0(?T,'^)). 
However, the aforementioned advantages of the Hermite and Laguerre-Fourier expansions are still 
valid and we leave the improvement of the conversion algorithm for the future work. 
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For both of the two approaches for deconvolution, we need some objective measurement scheme 
to assess the deblurred data in order to pick the best regularization parameter that gives the best 
deblurred image. With various values of the regularization parameter, we have a set of candidates 
for the estimation of the original image. In our work, we choose one with naked eye, which is not too 
blurry or too sharp. For a more systematic and objective assessment, a special kind of measuring 
method should be proposed. Furthermore, if we don't know the variance of the Gaussian motion 
distribution, which is a more practical situation, our tunable parameters are two: variance and 
regularization parameter. We also leave the development of an objective tuning algorithm for the 
future. 

6 Appendix 

On a polar plane and a circle, the Gaussian distribution densities are 




and 




k= 



oo 



Similarly we define the functions in SE(2) as 




and 




oo 



We want to show that 



Proof 




SE(2) 
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and 



(F2 * F,){g) = / F2{h)F,{h-' o g)dh. 

7SE(2) 



By changing the variable k — h ^ o g, we have 

(Fi * F2)ig) = [ F,{go k-^)F2{k)dk = j F2{h)F,{g o h-^)dh. 

JSE(2) JSE(2) 



g and h can be parameterized as 



cos ^ — sin ^ r cos 
sin 9 cos ^ r sin 
1 



and /i = h{R, O) 



cos © — sin © R cos $ 
sin © cos © i? sin $ 
1 



The multiphcation is 

h'^ og : 



goh ^ 



cos{e - 0) - sin(^ - 0) r cos(0 - Q) - R cos($ - 0) 
sin(^-0) cos(^-0) rsin(0-0) -i?sin($-0)) 
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cos(6l - 0) -sin(6'-0) r cos0 - i?cos(6' + $ - 0) 
sin(6' - 0) cos(6' - 0) r sin - sin(6' + $ - 0) 
1 



Therefore, 



{Fi*F2){g) = / F2{h)F,{goh-')dh 
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SE(2) i STT^f;:-^-' 27ri? 

Integration over © gives 
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Using the fact that the 5{R)/{2'kR) is a special delta function on a polar coordinate at singularity 
( = ), we have 



A;=— oo 



Similarly, 



SE(2) 



l-E 



(F2 * F^){g) = j F^{h)F^{h-' o g)dh 



e 



27ri? 



^ 00 



k=—oo 



Therefore, we showed that 
where 



k=—oo 



g = gir,(f),9) 



cos 9 — sin^^ rcos0 
sin 6 cos 6 r sin (p 
1 
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